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Abstract 

The study of biochemical pathways usually focuses on a small section of a protein 
interactions network. Two distinct sources contribute to the noise in such a system: 
intrinsic noise, inherent in the studied reactions, and extrinsic noise generated in 
other parts of the network or in the environment. We study the effect of extrin- 
sic noise entering the system through a nonlinear uptake reaction which acts as a 
nonlinear filter. Varying input noise intensity varies the mean of the noise after the 
passage through the filter, which changes the stability properties of the system. The 
steady-state displacement due to small noise is independent on the kinetics of the 
system but it only depends on the nonlinearity of the input function. 

For monotonically increasing and concave input functions such as the Michaelis- 
Menten uptake rate, we give a simple argument based on the small-noise expansion, 
which enables qualitative predictions of the steady-state displacement only by in- 
spection of experimental data: when weak and rapid noise enters the system through 
a Michaelis-Menten reaction, then the graph of the system's steady states vs. the 
mean of the input signal always shifts to the right as noise intensity increases. 

We test the predictions on two models of lac operon, where TMG/lactose uptake 
is driven by a Michaelis-Menten enzymatic process. We show that as a consequence 
of the steady state displacement due to fluctuations in extracellular TMG/lactose 
concentration the lac switch responds in an asymmetric manner: as noise inten- 
sity increases, switching off lactose metabolism becomes easier and switching it on 
becomes more difficult. 

Key words: Genetic switches and networks, noise in biological systems, nonlinear 
input function 
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Fig. 1. Steady-state displacement due to noise passing through a nonlinear uptake 
reaction. Extrinsic noise Pt enters the system through the uptake reaction which 
acts as a nonlinear filter hQ: for example, when the input signal has a Gaussian 
distribution, the output signal would have a skewed distribution. For small noise, 
varying the intensity of the input noise, without changing its mean P, varies the 
mean of the output noise (h(P t )). This causes a shift A(P) of the steady states 
of the system. The displacement direction depends on the shape of the filtering 
function h{P t ). The small-noise expansion method enables qualitative prediction of 
the displacement direction, by inspection of the following experimental data (here, 
an example sketch): the graph of the system's steady states vs. P and the graph of 



the u ptake rate vs. P. Real experimental data can be found e.g. in (|Qzbudak et al 
20o3 ). 



1 Introduction 



When studying biochemical pathways, we usually focus on a very small section 
extracted from a large network of protein interactions. Two distinct sources 
contribute to the noise in the studied system: Intrinsic noise due to the ran- 
domness of molecular collisions which result in chemical reactions of interest, 
and extrinsic noise, generate d in other parts of the network or in the external 
environment (iShibatal . 120051 ) . 



Computati onal methods developed based on the Gillespie algorithm (IGillespie 
19761 . 120071 ) enormously increased the popularity of theoretical studie s of in- 



trinsic s tochasticit v in biochemical reactions (see the most cited papers: iThattai and Oudenaarden 
f l200lh : iRao et al.1 (|2002h : lOzhudak et~ai1 rt2002T i). On the other hand grow- 



ing attention is being paid to the study of extrinsic noise (ISwain et al 



Rai and van Oudenaardenl . 120081 : IShahrezaei and Swainl . 120081 : IShahrezaei et al. 



2002 



20081 : iLei et all 120091 ) . in particular to the problem of discrimination between 



the effects of intrinsic and extrinsic noise components in studied systems 



2 



( Elo witz et al.l.l2002l; IRaser and O'Sheal . 12004 ; iPaulssonl . 12004 ; iPedraza and Oudenaarden 
200,4 iNewman et all l2006) . Experiments on lac gene in E. coli (jElowitz et al 



2002) and PH05 and GAL1 genes in Saccharomyces cerevisiae (IRaser and O'Shea 
2QQM) have demonstrated that the contribution of extrinsic noise to gene ex- 



pression may be significant. 

Noise propagates across a gene network by passing from one subsystem to the 
other through reactions which link particular sections of the network, or which 
link the network with the external environment. In this paper, we study the 
effect of extrinsic fluctuations entering the studied system through a nonlinear 
uptake reaction which acts as a nonlinear filter. In this way, one of the reaction 
rates in the system depends i n a nonlinear way on a fluctuating extrinsic vari- 
able. IShahrezaei et al.l ( 120081 ) simulated (with a Gillespie algorithm modified 
by inclusion of a randomly varying reaction rate) a simple general model of 
gene expression where one of the reaction rates was an exponential function 
of a slowly fluctuating Gaussian process. The resulting reaction rate had an 
asymmetric, log-normal distribution. This particular shape of the nonlinear fil- 
ter was chosen to conform the experime n tally measured di s tribut ions of gene 
expression rates (IRosenfeld et all 120051 ) . IShahrezaei et al.l ( 120081 ) found that 
the presence of this kind of noise shifts the mean concentration of the reaction 
product, however they did not explain that effect with analytical calculations. 



While IShahrezaei et al.l (120081 ) studied numerically the effect of extrinsic fluc- 
tuations whose average lifetime was longer than the shortest timescale of the 
system (thus interfering the relaxation dynamics of the system's variables), we 
focus on the analytical study of the effect of weak and rapid extrinsic noise, 
whose instantaneous fluctuations are too fast to directly interfere the system's 
dynamics. The assumption of wea kness and rapidi ty of the noise allows to use 
the small-noise expansion method (jGardineii . ll983l ). A measure for the rapidity 
of th e random fluctuations is the correlation time ( lHorsthemke and Lefeverl . 
1984 ). When the characteristic time scale of the noise (the correlation time) is 
much faster than the time scale of the reactions within the system of interest, 
then the noisy inpu t signal contributes to the kine tics of the system as an 
average signal only (lHorsthemke and Lefeverl . Il984 ). One could expect that, 
in such a case, varying the input noise intensity does not influence the be- 
havior of the system. However, the non-linear filtering function can transform 
the output noise distribution in such a way that its mean varies as the input 
noise intensity is varied (see Fig. []]). The system may therefore respond to 
extri nsic noise by shift i ng its steady state s with respect to the d eterministic 
ones dOchah-Marcinekl . 120081 : IRoccoI . 120091 : lOerstung et~all 120091 ) . 



A parti cularly interesting case is a multistable system, for example a bistable 
switch (lLaurent and Kellershohnl . ll999l : lFerrelll . l2002l ; IWilhelm . l2009h . a chem- 
ical reaction system where two distinct steady states are possible under the 
same external conditions. In the presence of a noisy input signal, the bistable 
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region can move to a concentration range where otherwise only one steady 
state was present, or bistabili ty can disappear for co ncentrations at which the 
system was initially bistable ( Ochab-Marcinek . 20081 ). Thus, as noise intensity 
increases, the existing steady states of the bistable system can be stabilized or 
destabilized, which results in the increase or decreas e of the escape time from 
one steady state to the other (jOchab-Marcinekl . 120081 ) . Th e switching time in a 
bistable genetic switch has been recently investigated by ICheng et al.l (120081 ) , 
as the measure for the robustness of the cellular memory to intrinsic noise. The 
simplicity of the system (one-dimensional equation of kinetics) allowed them 
for analytical calculation of the switching t ime as the mean first-passage time 
deriv ed from the Fokker -Planck equation (jGardinerl . 1 19831 ) . At a fixed noise 
level. ICheng et al.l (120081 ) studied the de pendence o f the switching time on the 
transcription rate. On the other hand, iLei et al.l ( 120091 ) simulated a similar 
system with a transcription delay. Extrinsic noise was the re applied in the 



same log-normal form as in the above-mentioned paper by IShahrezaei et al. 



(200 81) and the noise was strong enough to allow for bi-directional switching. 
Lei et al.l (120091 ) showed that the switching frequency decreases as the delay 
in the negative feedback loop increases. 



The study of stochastic control of metabolic pathways from the perspective 
of the small-noise expa nsion is relatively novel. It was initially presented in 
(jOchab-Marcinekl . 120081 ) where we used the small-noise expansion to calculate 
the noise-induced steady-state displacement for an arbi trary filtering function 
and t ested it on the reduced Yildirim-Mackey model (lYildirim and Mackeyl . 
20031 ) of the lac operon. Using a numerical simulation of that model, we showed 
that the induction time of the operon increases and the uninduction t ime de- 
creas es as a consequence of the noise-induced displacement. Recently, (jRoccol . 
20091 ) studied fluctuations in enzyme activity in Michaelis-Menten kinetics, 
expressing the noise-induce d steady-state d i splace ment as the Stratonovich 



drift. On the other hand, (IGerstung et all 120091 ) used the small-noise ex 



pansion method to study the steady-state displacement of concentrations of 
cis-regulatory promoters due to intrinsic fluctuations in transcription factor 
concentration combined with extrinsic fluctuations in the transcription factor 
synthesis rate. The non-linear input function (the filter) was the probability 
of promoter occupation in the form of h{n) = nj [K + as a function of the 
fluctuating number n of free transcription factors. 



While in (jGerstung et all 120091 ) the above form of the filtering function arised 
from the stationary solution of Master equation for the probability of finding 
n free transcription factors, in my present work the filtering function of the 
same form is the uptake rate of a certain substance P (the input signal), which 
enters the system through a catalytic Michaelis-Menten reaction (see Fig. H]). 
In a reaction that follows Michaelis-Menten kinetics, the reactant combines 
reversibly with an enzyme to form a complex. The complex then dissociates 
into the product and the free enzyme. If the total enzyme concentration does 
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not change over time, and the concentration of the substrate-bound enzyme 
changes much more slowly than those of the product and substrate (the quasi- 
steady-state assumption), then the intermediate enzymatic steps can be re- 
duced in the equations of chemical kinetics. The rate of product formation 
(P uptake) is then proportional to h(P) = P/(Km + P), where P is the 
concentration of P (lAtkinsl , 1 19981 ) . 



We show that, when the P concentration P* fluctuates with a Gaussian distri- 



bution , the n h{Pt) ha s an as ymmetric distribution, similarly as in lShahrezaei et al 
( 2008 ) and Lei et al. fj2009h . A simple calculus argument, based on the small- 
noise expansion method, allows to show that the noise-induced steady-state 
displacement depends solely on the shape of the nonlinear input function 
( Gerstung et al. . 20091 ). Consequently, we argue that the Michaelis-Menten 
type uptake function (monotonically increasing and concave) always generates 
the same type of the steady-state displacement. We show that this argument 
enables qualitative predictions of the steady-state displacement only by in- 
spection of experimental data, i.e. a) the graph of the uptake function vs. the 
mean of the input signal and b) the graph of system's steady states vs. the 
mean of the input signal. 



We show how, as a consequence of the steady-state displacement, the induc- 
tion/uninduction time of the lac operon changes depending on the intensity of 
TMG/lactose fluctuations. To illustrate the universality of t he phenomenon 



we co mpare the results for the Ozbudak model of lac operon (IQzbudak et al. 



20041 ). where the extracellular TMG uptake rate is an experimentally fitted 
function ~ P os (increa sing and concave ) , with the results for the reduce d 
Yildirim-Mackey model (lOchab-Marcinekl . 120081 : lYildirim and Mackevl . l200ah , 
where the extracellular lactose uptake rate is given by the Michaelis-Menten 
formula Pj (Km + P)- For both models we show how the Gaussian distribution 
of the input signal generates the skewed distribution of the Michaelis-Menten 
uptake rate. We analyze the validity range of the small-noise expansion in both 
cases. We compare the noise-induced steady state displacement and the noise- 
induced changes in the induction/uniduction time for the Ozbudak model 
with the previous l y pub lished results for the reduced Yildirim-Mackey model 
(jOchab-Marcinekl . l200sh . 



2 Theory 



2.1 Small-noise expansion method 



The input concentration of a substance P is modeled by the stochastic process 
Pt with the mean P and the variance a 2 . The passage through the uptake 
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reaction generates the output noise h(P t ). 

The equations of kinetics of the studied system are: 

£ = F(X,h(P t )), (1) 

where X is the vector of concentrations of substrates, and the system depends 
on P t through the function h(Pt) only. 



Assume that the noise is weak and slower than the time scale of the filter 
response but faster than the characteristic time scale of the system jTJ). Then 
the system experience s the output noise (after the pas sage through the filter) 
as its mean, (h(P t )) ( Horsthemke and Lefever . 19841 ). If the mean value of 



the output noise differs from the deterministic value of h(P), then the steady 
states of 

X = F(X,{h(P t ))) (2) 
are shifted with respect to steady states of the deterministic system 

X = F(X,h(P)). (3) 

It is possible to approximately find the displacement without the knowledge of 
the equations of kinetics, just by transformation of the X vs. P graph. While 
the deterministic system has stationary states X*, the stochastic system has 
quasi-steady states {X), around which its trajectories fluctuate (assuming that 
a noise-induced transition between multiple steady states is very unlikely). 
In the stochastic system, the steady states differ from the deterministic ones. 
This difference can be graphically shown as the shift of the graph of the steady 
states X* vs. P along the P axis0 (see Fig. by a function A(P), such that 

(X(P t ))=X*{P + A(P)). (4) 

The system depends on P t only in the function of the output process h(P t ), 
so 

(X(h(P t )))=X*(h(P + A(P))). (5) 
As assumed above, the system only experiences the mean of the output process 
(Fig. ED: 

(X(h(P t )))=X*((h(P t ))), (6) 



1 There are two reasons to consider the horizontal shift (with respect to P) and 
not the vertical one (with respect to (X)): a) The input signal Pt is the primary 
source of noise so it is easier to perform the calculations based on the dependence of 
the steady states on the variations in P t than on the variations of X, which anyway 
depend those of Pt. b) We will analyze the noise-induced shift in bistable systems, 
which can have two steady states for one value of P, so the vertical shift would be 
ambiguous. 
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And therefore, from the combination of Eqs. ((5j) and (jSJ): 



(h(P t )) = h(P + A(P)). (7) 

In other words, when the input process P t is stochastic (fluctuating around 
the mean P), then the mean of the output process h(P t ) is equal to the de- 
terministic output under the shift A(P) of the input flux. 

Assume that A(P) depends on P so weakly that it can be treated as a constant 
A. Expand the output process h(P+A) around the mean of the input process: 

(h(P t )) « h(P + A) = h(P) + h'(P)A + ... (8) 



(As shown in the Results section and in ( Qchab-Marcinekl . 2008 ). this approx 



imation works well for the example systems under study. In case when A(P) 
depends steeply on P, the approximation may break down.) 

On the other hand, expansion of (h(P t )) around P with respect to a fluctuation 
SP of the input process (P t = P + 5P at a given time t) yields: 

(h(P t )) = (h(P + 6P)) = h(P) + ti(P)(5P) + h"(P)(5P 2 ) + ... (9) 

But the mean deviation from the mean of the input process (SP) = and 
the variance of the input process (SP 2 ) = a 2 . Combination of Eqs. jEJ) and 
([9j) yields the approximate formula for the noise-induced shift of steady states 
with respect to deterministic steady states: 

h"(P) o , s 

A < p > = mrf m 



2.2 Michaelis-Menten uptake rate as the non-linear filter 



In the case of a s ystem where th e filtering function is generated by Michaelis- 



Menten kinetics (lAtkinsl . [1998), the small-noise expansion method is valid 



when the noise is slower than the intermediate enzymatic reactions, which had 
been reduced into h(P). Then, the enzymatic reactions are not perturbed by 
the fluctuations in P concentration, and the concentrations of the reactants can 
reach their steady state while the P concentration is approximately constant. 



The noise correction A(P) to the Michaelis-Menten rate h(P) = P/(K M + P) 
can be viewed as a correction to the Michaelis constant Km- The displacement 
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2 4 6 8 10 
Extracellular sugar concentration P [mM] 

Fig. 2. The TMG uptake rate (fTT|) for the Ozbudak model (solid line) and the lactose 
uptake rate (fl~9l) for the Yildirim-Mackey model (dashed line). The mathematical 
formulas for the uptake rates are different in both models, but as they describe the 
same uptake mechanism, their graphs have a similar Michaelis-Menten shape with 
hi > and h" < 0, which generates the same type of the steady-state shift due to 
noise in extracellular TMG/lactose concentration (see the Results section). 

P + A(JP) = P(l + A(P)/P) is equivalent to the correction 



K 



M,noise 



K 



M 



K 



M 



1 + 



MP) 
p 



1 - 



(11) 



(K+P)P 



When the input signal P t has a Gaussian distribution with the mean P and 
variance a 2 , the output signal h(P t ) = P t /(K+P t ) (Fig. has an asymmetric 
probability density function (Fig. [3]b): 



p(h) 



K 



V2 



exp 



Kh 



2a 2 \ l-h 



-P 



:i2) 



(See the Appendix [A] for the details of the calculation.) The mean (h(P t )) = 
J hp(h)dh/Af of (Tl2l) exists if > h > 1, which corresponds to > P > +oo. 
The normalization constant J\f = (l/2)(erf [Py2/(2a)] + 1), and the accurate 
values of (h(P t )) can be computed numerically. 



3 Models 



The lac operon is one of the most extensively studied examples of a bistable 
genetic switch. In a certain range of extracellular TMG/lactose concentration, 
the gene transcription can be in one of two discrete states, either fully in- 
duced or uninduced. Bistability is generated by the positive feedback loop, 
where high lactose/TMG concentration in the cell causes weak repression of 
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Fig. 3. Probability density functions for the TMG/lactose uptake rates when the 
fluctuating concentrations of the extracellular TMG/lactose have a Gaussian dis- 
tribution with the variance a 2 . Vertical lines indicate the mean, which is slightly 
shifted with respect to the deterministic uptake rate (solid line), a) Ozbudak model, 
b) Yildirim-Mackey model. 

the lac gene transcription and, in consequence, more permease is produced 



(Huber et al.. 


1980; 


Wrieht et al.. 


1981: 



1991 



Ozbudak et all l2004h the rate of the TMG/lactose uptake into the cell 



is of Michaelis-Menten type. Extracellular lactose concentration experienced 
by E. coli living in its nat ural environment may fluctuate, due to high mobility 
of the bacterium (see e.g. IPiLuzio et al.l ( 120051 )). granularity of the intestinal 
content and motions of intestinal villi. 



3.1 Ozbudak model 



The model ( Ozbudak et al. . 20041 ) is based on the experiment where a gene 
coding for a fluorescent protein was incorporated under the control of the lac 
promoter in E. coli bacteria, in order to indicate the transcription activity 
of the lac gene. The equations of chemical kinetics for this system are the 
following: 

R - 1 

ifr ~ 1 + (x/x ) 2 (13) 
dy 1 

T «Tt =a TTRlR - y (14) 

(1 T 

r x — = p G h(T)y - x (15) 

x denotes the intracellular TMG concentration, y is the concentration of per- 
mease in green fluorescence units, is the total concentration of the repres- 
sor, and R is the concentration of active repressor. The active fraction of the 
repressor is a function of the TMG concentration x, with half-saturation con- 
centration Xq. a is the maximum rate of generation of permease, Rq is the 
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half-saturation of the repressor. Permease is depleted in a time scale r y , due 
to a combination of degradation and dilution. TMG enters the cell at the rate 
h(T) (Fig. [2|) proportional to y and to the glucose uptake rate /3q, and is 
depleted with time constant r x . 

Steady states of y (in xq units, further on called the 'scaled units') are given 
by: 

i + (p G h(T) y y 

l + f + (/W>) 2 
T is measured in /iM. In a certain range of TMG concentrations, the cubic 
equation (fT6l) has three roots (two stable fixed points separated by one un- 
stable fixed point), which very well reproduces the experimentally observed 
switch-like behavior of lac operon. In this study, we assume that the ex- 
tracellular glucose concentration is zero. Then, the system is bistable for 
3.39 < T < 24.40 /iM. We have chosen the value of r x — 1 min significantly 
larger than the time scale of the noise ( 1201) , and at the same ti me much smaller 



than r y = 216 min, which conforms the experimental results ( IMettetal et al 



20061 ) reporting r x less than the measurement resolution of 35 min. 

The TMG uptake rate function 

h(T) = 1.23 x 10" 3 T ' 6 , (17) 



as well as the value s of the parame t ers (T able IC.ll) , have been fitted from the 



experimental data ([Ozbudak et all 12004 ). Since the formula (fT7]l is the result 



of fitting, it does not have the classical Michaelis-Menten form T/(K + T), 
but its graph has the characteristic Michaelis-Menten shape, increasing and 
concave. 



3.2 Yildirim-Mackey model 



The Yildirim-Mackey model (jOchab-Marcinekl . 120081 ; lYildirim and Mackeyi . 



2003h consists of three equations of chemical kinetics for mRNA (M), allolac 



tose (A) and lactose (L) concentrations in the E. coli cell: 

~W - aM l+K 2 R t ot+K x A* + 1 - IM M 



dA 
dt 



^M(a A ^ rL -(5 A ^)-~ lA A (18) 
f =k P M (a L h(L e ) - /3 L -a A k B M ^ - % L 
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Ozbudak model: 
Bistable region 



Y-M model: 

Bistable region 

-0.4 1 1 ' 1 1 ' 1 

10 20 30 40 50 60 70 

Extracellular lactose/TMG concentration [(iM] 

Fig. 4. Comparison of the steady-state displacements calculated using the small-noise 
expansion method for the Ozbudak model (a) and the Yildirim-Mackey model (b). 
Rectangles indicate bistable regions of both models. In the Ozbudak model, the 
noise-induced shift of steady states is weaker than in the Yildirim-Mackey model and 
it depends more strongly on the extracellular sugar concentration. The displacement 
is the strongest on the left boundary of the bistable region and the weakest on its 
right boundary. 



a and (3 denote the gain and loss rates for the reactions. K x is the equilibrium 
constant for the repressor-allolactose reaction. K 2 is the equilibrium constant 
for the operator-repressor reaction, and R to t is the total amount of the repres- 
sor. The 7 = 7 + fi are the coefficients for the terms representing decay of 
species due to chemical degradation (7) and dilution (/z). Even if allolactose is 
totally absent, on occasion repressor will transiently not be bound to the op- 
erator and RNA polymerase will initiate transcription. Although the mRNA 
production rate dM/dt would be then nonzero (a leakage transcription), it 
is necessary to add an empirical constant T t o the model to obtain a leak - 
age rate that agrees with experimental values ( Yildirim and Mackey . 2003T ) . 
The allolactose gain and loss rates and the lactose loss rate depend on the 
/3-galactosidase concentration (an anzyme breaking down lactose into allolac- 
tose), which is proportional (kg factor) to the mRNA concentration. Similarly, 
the lactose gain and loss rates depend on the permease concentration, propor- 
tional (kp factor) to the mRNA concentration. See Table IC.2I for the values 
of parameters. 



The system is bistable for 27.7 fj,M 
function 

h(L e ) 



< L e < 61.8 fjM. The lactose uptake rate 



(19) 



has the Michaelis-Menten form, to c o nform the data of lLolkema et al.l (119911 ) ; 
Huber et ali (Il980l ): IPage and Westl rtl984T ): IWright et all ( jl98ll ). 
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3.3 Fluctuations 



The fluctuations in TMG/lactose concentration have been modeled by Ornstein- 
Uhlenbeck noise. For the small-noise expansion to be valid, the correlation time 
of the noise 

tou = 1.2 s (20) 

has been chosen significantly larger than the time s cale of the enzymat ic 
TMG/lactose uptake by permease (0.1 s, according to Wright et al. ( 19861 )). 
but smaller than the fastest time scale of the system: r s 



1 min 

for the Ozbudak model and r sys m 10 _1 min for Yildirim-Mackey model 



SYS 



dYildirim and Mackevl . l2~0~0~i l. Ta king into account the bacterial motility (mean 
velocity of E. coli is ~ 30//m/s, (jDiLuzio et all 120051 )). the granularity of the 
intestinal content and motions of intestinal villi, one can suppose that the 
fluctuation rapidity assumed here is realistic (jOchab-Marcinekl . 120081 ) . 



The details of the numerical simulation of the stochastic process are presented 
in the Appendix iBl 



4 Results 



4-1 Asymmetric distribution of the Michaelis-Menten uptake rate 



For the Mackey model, the probability density function of the Michaelis- 
Menten uptake rate h is given by the formula (fT2]) . For the Ozbudak Model, 
the probability density function of the uptake rate is described by a different 
formula (see the Appendix [A] for the details of the calculation): 



p(h) 



3&s/3 v /2^ 



: exp 



2a 2 




T 



(21) 



with b = 1.23 x 10~ 3 , because h(T) is given by the fitted function ( TT71) which, 
however, has a similar shape to the Michaelis-Menten function. Consequently, 
the probability density functions for both models also have similar asymmetric 
shapes, and for small noise their mean values decrease as the noise intensity 
increases (Fig. [3]). For larger noise, the mean values begin to increase (compare 
Fig. C]). 
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Fig. 5. Comparison of the steady-state displacements (1231 I24|) calculated using 
the small-noise expansion method with values obtained from the simulation of 
the Qzbudak model (a), and corresponding results for the Yildirim-Mackey model 
( Ochab-Marcinek . 20081 ) (b). The standard deviations of the fluctuations were 
a = 1. 5 uM for the Qzbudak mo del and a = 10 fjM. for the Yildirim-Mackey 
model ( Yildirim and Mackey . 20031 ). The response of both models to fluctuations in 
extracellular TMG/lactose concentration passing through a Michaelis-Menten-type 
uptake reaction is very similar: The steady states shift to the right. The displacement 
is larger for lower TMG/lactose concentrations and smaller for higher concentrations. 



4-2 Calculation of the steady-state displacement using small-noise expansion 



The steady-state displacement (flOl) does not depend on the kinetics of the 
system. It only depends on the input noise intensity and on the shape of the 
uptake function: its monotonicity and concavity. Therefore, only knowing the 
graph of the steady states vs. P , one can predict the changes in its stability 
due to noise passing through the uptake reaction. 



For a monotonically increasing and concave uptake function, such as Michealis- 
Menten (Fig. EJ), h' > and h" < 0, so the steady states always shift to the 
right (in the direction of higher concentrations of P): 

A(P) < (22) 

The formulas for TMG/lactose uptake rates (fTTj) . ( fT9]l are different in both 
models, but as they describe the same uptake mechanism, their graphs have 
a very similar Michaelis-Menten shape (Fig. EJ) with h! > and h" < 0, which 
causes a steady-state shift to the right. In the Ozbudak model, 

A(T) = (23) 
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Fig. 6. The values of the steady-state displacement, calculated using the small-noise 
expansion method, compared with the values obtained from simulations, a) the 
Ozbudak model, b) the Yildirim-Mackey model. The noise intensities were same as 
in Fig. 03 



while in the Yildirim-Mackey model (lOchab-Marcinekl . 120081 ). 

1 



A(L e 



(K Le +L e 



-a 



(24) 



Within the bistable regions, the steady-state displacement due to noise in 
the Ozbudak model is smaller than in the Yildirim-Mackey model (Figs. SI 
[5]). A common feature of both models is a larger shift for low extracellular 
TMG/lactose concentrations and a smaller shift for high concentrations. 

In the Ozbudak model as well as in the Yildirim-Mackey model, the results ob- 
tained analytically by small-noise expansion fl23l l24l) were in a very good agree- 
ment with the simulation results (Figs. [5J [6]). The behavior of the Ozbudak 
model was very similar to the behavior of the Yildirim-Mackey model: In both 
cases the graph of the steady states vs. the mean extracellular TMG/lactose 
concentration shifted to the right. The shift was larger for lower TMG/lactose 
concentrations and smaller for higher concentrations. 



4-3 Range of validity of the small-noise expansion 



For both models (within the given choice of parameters), we compared the val- 
ues of A calculated using different methods for chosen constant values T or L e 
as noise intensity and time scale of the noise were changed (Fig. Ej). The values 
of A were: a) obtained from the simulation, b) calculated accurately (the mean 
(h(P t )), computed using the distributions (TT2l) or (l2~Ti) . was substituted into the 
equations of kinetics instead of h(P t )), and c) calculated using the small-noise 
expansion. The results are consistent with those expected: The expansion (JHJ) 
is valid when A C P, and indeed, in both models the approximation breaks 
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Fig. 7. The range of validity of the small-noise expansion: Comparison of the 
steady-state displacement calculated using different methods (simulation, accurate 
calculation, small-noise expansion) for constant values T or L e as noise intensity and 
time scale of the noise are changed. Top panels: Varying noise intensities a 2 , a) Ozbu- 
dak model, f = 3.6 fjM, t u = 1.2 s = 0.02 min. b) Mackey model, L e = 30 /xM, 
tou = 1.2 s = 0.02 min. Bottom panels: Varying noise time scales tou- c) Ozbudak 
model, f = 3.6 fjM, a 2 = 1 \iM 2 . d) Mackey model, L e = 30 /xM, a 2 = 50 /iM 2 
(the accurate value overlaps with the approximated one). 



down when A/T or A/L e are greater than the order of 10~ 2 . Moreover, the 
time scale of noise should be less than the time scale of the system, and for the 
Yildirim-Mackey model the approximation is valid for t Q u < 0.1 min while 
the fastest charac t eristic time for the left bifurcation point was T sys = 0.4 min 
( jOchab-Marcinekl . 120081 ). For the Ozbudak model, the approximation breaks 
down at tou ~ O.lmin, while T sys = lmin. The results of the accurate calcula- 
tion differ slightly from the small-noise approximation because of the Taylor 
expansion cut-off. There is also a systematic difference, increasing with noise 
intensity, between the simulation results and those calculated accurately. This 
difference is due to the reflecting boundary conditions used in the simulations, 
which add a contribution from the trajectories reflected at P t = 0. 
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Fig. 8. Increasing noise intensity inhibits induction and accelerates uninduction in 
the studied models of the lac operon. Mean induction/uninduction time was mea- 
sured in simulations of the Ozbudak model (a) with different noise intensities and 
different mean TMG concentrations f: 24.5 fM (A), 24.6 fM (B), 24.7 fM (C), 
3.41 /jM (D), 3.40 fM (E), 3.39 fiM (F). These results ar e comp ared with the re- 
sults for the Yildirim-Mackey model (b) ( Ochab-Marcinekl . 20081 ). where the mean 
extracellular lactose concentrations L e were: 62 fiM (A), 63 fM (B), 65 fM (C), 
27.9 fM (D), 27.8 fM (E), 27.7 fiM (F). 



4-4 Induction/uninduction time 



We compared the mean induction/uninduction times for both models, ob- 
tained from the simulations with different noise intensities and different TMG 
concentrations (Fig. [BJ). Both models are robust to fluctuations in extracel- 
lular TMG/lactose concentration. Noise-induced switching from the induced 
to uninduced state due to nois e driving the system out of the steady state 
(IHorsthemke and Lefeverl . 1 19841 ) was possible only for concentrations close to 
the boundaries of the bistable region. 



Therefore, to study the uninduction we performed simulations in which the 
system started within the bistable region, in the induced state very close to 
the left bifurcation point (Tab. IC.3I ). When the noise intensity is zero, the 
system remains in the steady state and the uninduction time is infinite. But 
as the noise intensity increases, the bistable region shifts in the direction of 
larger TMG concentrations, effectively destabilizing the induced state and 
thus decreasing the mean uninduction time (Fig. [9K). 

On the other hand, within the bistable region the noise-driven switching from 
the uninduced steady state to the induced one was impossible in the range 
of the noise intensities for which the small-noise expansion method is valid 
(Fig. [9b). To observe the transitions to the induced state, we had to set the 
initial conditions outside the bistable region (Tab. IG3I) . The unstable initial 
positions become closer to the bistable region, and thus become more stable, 
as this region shifts in the direction of larger TMG concentrations. The trajec- 
tories spend more time in the vicinity of the starting point before they switch 
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Fig. 9. Increasing noise intensity inhibits induction and accelerates uninduction in 
the Ozbudak model of the lac operon: Example trajectories for different initial states 
of the system: a) Induced, left bifurcation point (T = 3.39 ^M). b) Uninduced, right 
bifurcation point (T = 24.4 /j,M). c) Uninduced, initial concentrations same as the 
coordinates of the right bifurcation point, but the TMG concentration is slightly 
beyond the bistability region (T = 24.5 /iM). 

to the induced state, so the mean induction time increases as noise intensity 
incr eases (Fig. Eb) . T h is effe ct is same as in the Yildirim-Mackey model (Fig. 
ISlb. lOchah-MaTcinekl (l2008h ). 



5 Discussion 



In spite of the differences between the kinetics of the Ozbudak and Yildirim- 
Mackey models, and differently defined TMG/lactose uptake rate functions 
(the noise filters), the behavior of both models is qualitatively the same. The 
steady-state shift due to small noise and the consequent changes in the in- 
duction/uninduction times depend on the characteristics (monotonicity, con- 
cavity) of the filtering function and not on the details of the kinetics of the 
system. 

Gaussian fluctuations of the extrinsic input signal, which enter the system 
through the Michaelis-Menten uptake reaction, generate an asymmetric dis- 
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tribution of its rate. This effect is similar as in IShahrezaei et al.1 (120081 ) and 
Lei et al.1 (l2009h . where Gaussian no ise enters the syst e m thr ough an exponen- 



tial function. (Note that, however, in Shahrezaei et al. (120081 ) the noise was too 
slow (slower than the fastest time-scale of the system) and in lLei et al.l (120091 ) 
the noise was too strong for the small noise-expansion to be valid.) While in the 
latter cases the rate distribution is log-normal, the distribution the Michaelis- 
Menten uptake rate has an opposite skewness. Due to the asymmetry of the 
distribution, the mean uptake rate varies as the input noise intensity is varied. 
This gives rise to the shift of the steady-state concentrations of the studied 
reaction system. 



The small-noise expansion turns out to be a useful tool for prediction of the 
steady-state displacement due to weak and rapid noise. Within the validity 
range of the method, the distribution of the input noise is not important (it 
can be non-Gaussian as well): It is only the shape of the uptake function that 
determines the direction and size of the steady-state displacement for a given 
input noise intensity. The simplicity of this approximate calculation allows for 
qualitative predictions, based on experimental data only, of the response of 
the studied system to extrinsic noise passing the uptake reaction of a given 
type. Assuming that we can only measure the intensity of the input noise, all 
information needed to predict the steady-state displacement are: a) Experi- 
mental the graph of the system's steady states vs. the mean input signal and 
b) Experimental graph of the uptake rate function vs. the mean input signal. 
Even if this data is not precise or already includes noise (i.e. is shifted due to 
noise), it enables qualitative predictions of the steady-state displacement when 
noise intensity increases or decreases. The shape (monotonicity and concavity) 
of the uptake rate function indicates the direction of the noise-induced shift 
of the stationary states, whose approximate positions are known from the ex- 
perimental graph. In particular, when the uptake rate is of Michaelis-Menten 
type (monotonically increasing and concave), the graph of the system's steady 
states vs. mean input signal always shifts to the right. In case of a bistable 
genetic switch, one can then predict that, as a consequence of such a shift, the 
extrinsic noise can selectively facilitate or inhibit induction and uninduction. 

In principle, it would be possible to measure the shift experimentally on the 
level of the Michaelis-Menten reaction rate h(P) = P/(Km + P). The shift 
could be detected as the noise correction (fTTl) to the Michaelis constant K M) 
which could be read out from the experimental plot of h(P). However, in the 
studied models of the lac operon the noise would be probably too small for 
experimental detection of the shift of the reaction rate plot. But the computer 
simulations suggest that it would be easier to experimentally detect the shift 
by measurement of the switching times of the lac switch. 

The measurements of induction/uninduction times in the simulations of two 
example lac operon models demonstrate that the bistability of the lactose 
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utilization mechanism is robust to fluctuations in extracellular TMG/lactose 
concentration. However, even as small a noise as used in this study (to fulfill 
the conditions of validity of the small-noise expansion) can have a significant 
effect on the induction/uninduction time. The values by which it increases or 
decreases depend on the choice of model parameters and initial conditions. The 
effect is the strongest when the initial conditions are close to the bifurcation 
points (boundaries of the bistable region). For example in the Ozbudak model, 
when the initial state of the system is the left bifurcation point (induced state), 
then the uninduction time changes from infinity (at zero noise) to ~100 hours 
(when the standard deviation a of the TMG fluctuations is 1.5 /iM). On the 
other hand, when the initial state is the right bifurcation point (uninduced 
state), then the noise-driven induction is impossible. When the initial state is 
close to the right bifurcation point, but out of the bistable region, the induction 
time is increased by noise. For example, when the mean TMG concentration is 
24.5 fiM, then at zero noise the induction time is ~190 hours, but at the noise 
of o = 1.7/iM the induction time is ^25 hours. These effects ar e qualitatively 
same as in the Yildirim-Mackey model (lOchab-Marcinekl . 120081 ) . 



Thus, the fluctuations in extracellular TMG/lactose concentration facilitate 
the switching off of the TMG/lactose metabolism, but at the same time 
they prevent the metabolism from switching on. This effect will be valid 
for any model of lac operon, provided that the TMG/lactose uptake rate 
is of Michaelis-Menten type. This suggests that in the presence of noise the 
possibility of random switching on the metabolism of TMG/lactose is more 
strongly protected than the possibility of random switching it off. One can 
speculate whether this protection against accidental induction of metabolism 
is connected with preventing an unnecessary energetic effort. To answer this 
question, one should analyze the lactose utilization system in E. coli from the 
energetic point of view. 



In other systems, a different direction of the noise-induced steady state dis- 
placement is possible. For example, when extrinsic noise enters the system 
th rough the exponential f uncti on (monot i nically increasing and convex), as 



m 



Shahrezaei et al.1 (120081 ) and iLei et al.1 (120091 ) . then one can predict that 



for weak and rapid noise the graph of the steady states vs. mean noise in- 
tensity will shift in the opposite direction than that of the systems with the 
Michaelis-Menten uptake. On the other hand, if the uptake rate is a Hill func- 
tion P n /(K + P n ) then the formula (TTOl) for the steady-state displacement can 
change its sign for different concentrations of P and, under certain conditions, 
bistability may emerge due to noise in place of a graded response. Similar 
effects are possible for systems wi th non-monotonic input functions gen erated 



by incoherent feed-forward loops (IKaplan et all 120081 ; iKim et all 120081 ) 



The analysis presented applies to weak and rapid noise from one dominat- 
ing external source. But even in the systems where the contribution of other 
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noises is present, the method may be of use to interpret the experimental mea- 
surements in terms of the discrimination between the effects of noises which 
originate from different sources. 
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A Transformation of probability density functions 



When h(P) is an either monotonically increasing or decreasing function of a 
random variable P, and P is given by the probability density function q(P), 
then the probability density function of h is given by the formula: 



p(h) = q(P(h)) 



dP{h) 



dh 



(A.l) 



where P(h) is the inverse function of h(P). The formula is obtained by the 
change of variables in the integration (see e.g. iMiller and Milien fl2004h ): 



Prob(/i(Pi) <h< h(P 2 )) = Prob(Pi < P < P 2 



(A.2) 



P> 



q{P)dP 



h(P 2 ) 
h(Pi) 



q{P{h)) 



dP{h) 



dh 



dh 



h(P 2 ) 
HPi) 



p(h)dh 



(The absolute value is needed when h(P) is decreasing. 



B Simulation details 



The fluctuations in the extracellular TMG/lactose concentration are mod- 



elled by the Ornstein-Uhlenbeck (OU) process (IGardineii . ll983l ) with reflecting 
boundary conditions in P t = 0. £(£) is a Gaussian white noise of intensity 7 
and autocorrelation (£(t)£(s)) = 5(t — s): 



dP 

-i = -e(p t -p) + rt(t), p>o. 



(B.l) 



The correlation time of the noise tou = 1/0. The variance of the OU process is 
a 2 = r y 2 /29. We assume a small noise intensity and P sufficiently far from the 
reflecting boundary, so that the contribution of reflected 'tail' of the Gaussian 
distribution can be neglected and we can use formulas for the unbounded OU 
noise for the mean and variance. The noise intensity is varied in the simulations 
by varying the value of 7. . Numerical in tegration of t he equations of k inetic s 



has been done using the Euler scheme (IPress et all 119931 ; iMannellal . 120021 ). 
The timestep 5t = 2 ■ 10~ 3 min has been chosen significantly smaller than the 
time scales of the studied systems. 

To estimate the mean induction/uninduction time in the simulations, the time 
was measured until trajectories X(t) got into a close neighborhood of the 
other deterministic stationary state (of a radius D = \Jj2iSX~f = 0.1 [scaled 
units] for the Ozbudak model and D = 5 fj,M for the Yildirim-Mackey model 
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(lOchab-Marcinekl . 120081 )). The initial points of the trajectories X were the 
deterministic steady states within the bistable region, or points outside the 
bistable region, but close to its boundaries (see Tables |CT3~1 IC.4I) . The number 
of simulation runs for calculating the mean induction/uninduction time was 
N = 100 for the Ozbudak model and N = 1000 for the Yildirim-Mackey 
( lYildirim and Mackeyl . 120031 ) model. 



C Tables 



p = l 


i Rt 
^ Ro 


167.1 a 


a 




100.5 a 




= 0) 


100 a 


r y 




216 min b 






1 min c 



Table C.l 

Parameters of the Ozbudak model: a ) (lOzhnriak et all 120041 ). b ) dMettetal et al 



2006) , c ) chosen for this study within the range of 0..35 min reported by 
(jMettetal et~al1 . l200fih . 



r 


7.25 x 10" 7 mM/min 


/i 


0.0226 min" 1 


a A 


1.76 x 10 4 min" 1 


t b 


2.0 min 


a B 


1.66 x 10~ 2 min -1 


tm 


0.1 min 


a L 


2.88 x 10 3 min" 1 


TP 


0.83 min 


% 


9.97 x 10~ 4 mM/min 


K 


7.2 x 10 3 


ap 


10.0 min" 1 


K x 


2.52 x 10 4 mM" 2 


(3a 


2.15 x 10 4 min"l 


K A 


1.95 mM 


Pl 


2.65 x 10 3 min" 1 


K L 


0.97 mM 


7a 


0.52 min" 1 


K Le 


0.26 mM 


IB 


8.33 x 10~ 4 min" 1 


K Ll 


1.81 mM 


1L 


0.0 min" 1 


k B 


0.677 


1M 


0.411 min" 1 


kp 


13.94 


IP 


0.65 min" 1 







Table C.2 

Parameters of the Yildirim-Mackey model ( Yildirim and Mackey . 20031 ). 
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T[uM) 


X; [scaled units! 


Vi [scaled units! 


x /-[scaled units! 

J I J 


iif [scaled units! 


A: 24.5 a 


1.012 


1.210 


82.23 


98.09 


B: 24.6 a 


1.012 


1.210 


82.44 


98.10 


C: 24.7 a 


1.012 


1.210 


82.65 


98.12 


D: 3.39 


13.228 


51.700 


0.1577 


0.6163 


E: 3.40 


13.707 


53.475 


0.1580 


0.6164 


F: 3.41 


14.012 


54.569 


0.1583 


0.6164 



Table C.3 

Initial (i subscript) and final (/ subscript) concentrations of intracellular TMG (x) 
and permease (y) for the simulation measurements of mean switching time in the 
Ozbudak model. Trajectories A, B, C started from outside the bistable region, but 
very close to its boundaries. Initial points for these trajectories are the coordinates 
of the right ( a ) bifurcation point.) 



L e \jM\ 


Ai\fM] 




Mi[10- 2 /iM] 


A f [»M] 


L f \pM] 


M f [lQ- 2 ^M] 


A : 62.0 a 


14.1 


224 


0.360 


393 


285 


80.8 


B : 63.0 a 


14.1 


224 


0.360 


399 


289 


82.5 


C : 65.0 a 


14.1 


224 


0.360 


412 


298 


85.9 


D : 27.9 


109.0 


129 


9.406 


4.00 


94.6 


0.212 


E : 27.8 


104.0 


129 


8.600 


4.00 


94.0 


0.212 


F : 27.7 


96.9 


128 


7.521 


3.98 


93.7 


0.212 



Table C.4 

Initial (i subscript) and final (/ subscript) concentrations of allolactose {A), lactose 
(L), and mRNA (M) for the simulation measurement s of mean switching time in the 
Yildirim-Mackey model ( Yildirim and Mackey . 20031 ). a ) Trajectories starting from 
outside the bistable region, but very close to its boundaries. Initial points for these 
trajectories are the coordinates of the right bifurcation point. 
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